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The motion of one-hundred point vortices in a circular cylinder is simulated numeri- 
cally and compared with theoretical predictions based on statistical mechanics. The novel 
aspect considered here is that the vortices have greatly different circulation strengths. 
Specifically, there are four strong vortices and ninety-six weak vortices, the net circula- 
tion in either group is zero, and the strong circulations are five times larger than the weak 
circulations. As envisaged by Onsager [Nuovo Cimento 6 (suppl.), 279 (1949)], such an 
arrangement leads to a substantial amplification of statistical trends such as the preferred 
clustering of the strong vortices in either same-signed or oppositely-signed pairs, depend- 
ing on the overall energy level. To prepare the ground, this behaviour is illustrated here 
first by a simple toy model with exactly solvable statistics. A microcanonical ensemble 
based on the conserved total energy E and angular momentum M for the whole vortex 
system is then used, in which the few strong vortices are treated as a subsystem in con- 
tact with a reservoir composed of the many weak vortices. It is shown that allowing for 
the finite size of this reservoir is essential in order to predict the statistics of the strong 
vortices accurately. Notably, this goes beyond the standard canonical ensemble with pos- 
itive or negative temperature. A certain approximation is then shown to allow a single 
random sample of uniformly distributed vortex configurations to be used to predict the 
strong vortex statistics for all possible values of E and M. Detailed predictions for the 
energy, two-vortex and radial distribution functions of the strong vortices are then made 
for comparison with three simulated cases of near-zero M and low, neutral, or high E. 
It is found that the statistical mechanics predictions compare remarkably well with the 
numerical results, including a prediction of vortex accumulation at the cylinder wall for 
low values of E. 
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I. INTRODUCTION 

The application of statistical mechanics to two-dimensional point vortex dynamics 
was first suggested by Onsager^ in a landmark paper in 1949, in which he sketched a 
possible explanation for the formation of coherent vortices on statistical grounds linked 
to the possibility of negative (statistical) temperatures for point vortex systems in a 
bounded domain. Onsager's suggestions have continued to attract some interest because 
two-dimensional fluid dynamics is a relevant paradigm in many applications such as 
geophysical fluid dynamics (in which the large-scale quasi-horizontal flow along stratifi- 
cation surfaces in the atmosphere or oceans is approximately layerwise two-dimensional), 
or plasma dynamics imder certain conditions^. Indeed, a successful application of statis- 
tical mechanics to problems in these fields suggests a hidden degree of predictability that 
could easily be obscured by conventional direct numerical simulations. Such statistical 
predictability can be exploited quantitatively at much lower computational cost than 
by the brute force simulations of ensembles of many individual flow realizations. This is 
a compelling vision if one considers, for example, the hugely expensive geophysical cli- 
mate and weather simulations, for which currently only a very small ensemble of direct 
simulations at very low resolutions is feasible. Statistical theories based on Onsager's 
ideas (and on others that go beyond the point vortex idealization) have already been 
used successfully for predicting the detailed behaviour of certain idealized geophysical 
flow problems^. Formidable obstacles remain in order to make such theories applicable 
in practice, but the potential rewards are great. 

Onsager described his ideas only in qualitative form and the detailed theoretical ex- 
ploration of these issues only began with the development of a mean-field theory by 
Montgomery and Joyce*'® which has since been extended and refined mathematically in 
many ways over the years^"^. Accurate direct numerical simulations of point vortices 
over long times have become feasible over the last two decades^°'^^. It was suggested^^ 
that for geophysical applications the most relevant values for the number of vortices N 
lie in an intermediate range between the regime of low-dimensional chaos (where N may 
be less than 10) and the 'thermodynamic' regime in which iV ^ oo in some subtle limit. 
This point of view is also taken here, where only cases with N = 100 are studied. 

Now, the topic of the present paper goes back to a suggestive remark made by Onsager 
concerning the characteristic appearance of vortex distributions in a negative tempera- 
ture state, and which to my knowledge has not been considered explicitly since. In his 
masterful succinct exposition, Onsager says that in such a state 

"...vortices of the same sign will tend to cluster, — preferably the strongest ones — , 
so as to use up excess energy at the least possible cost in term of degrees of freedom 
. . . the weaker vortices, free to roam practically at random, will yield rather erratic and 
disorganised contributions to the flow" (my italics). 

This encapsulates two crucial insights: (a) that strong vortices will be more predictable 
than weak ones; and (b) that the maximum disorder of the flow as a whole (which 
is implied by the statistical theory) will be achieved in a remarkably inhomogeneous, 
composite manner, in which the relative order of the strong vortices will be more than 
compensated for by the increased disorder of the weak ones. The present paper is an 
attempt to verify and exploit both of these insights in the simplest possible setting that 
allows comprehensive numerical and theoretical exploration. 

To the best of my knowledge, detailed previous studies have all focused on the case 
of identical (or nearly identical^ "'^) absolute vortex circulations, or on the even more 
restricted case of identical vortex circulations throughout. In these cases, the statistical 
behaviour of the flow is in some sense directly determined by the constraint of fixed total 
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energy: large energy values must require the coming together of same-signed vortices, 
and vice versa for low energy values (see ^ below). In other words, the occurrence of 
vortex clusters is then enforced directly by the total energy constraint. However, once 
different vortex strengths are present there is scope for interestingly different behaviour, 
in which some vortices cluster spontaneously whilst others do not. 

The plan for the paper is as follows: §|| introduces the main features of the studied 



Hamiltonian point vortex system and corrects some minor errors in related works; §11] 
briefly illustrates Onsager's suggestion by an analogy with a simple toy model, which 
serves to prepare to ground for the vortex theory; §^ presents direct numerical simu- 
lations of the vortices; derives statistical mechanics predictions and applies these to 
the simulations; and some concluding remarks are given in §VI. 



II. HAMILTONIAN POINT VORTEX DYNAMICS 

The Hamiltonian form of the equations of motion for N point vortices with circulations 
Ti and instantaneous Cartesian coordinates Xi{t) — {xi, yi) (where i = 1, 2, . . . , N) are: 

r,^ = +^, r,^ = -— . (ii.i) 

At dyi ' At dxi 

The pairs {xi, yi) are canonical phase space coordinates, with invariant phase space vol- 
ume element da;ida;2 • • -Axm- In the special case of a circular cylinder with radius R 
centred at the coordinate origin^ the wall boundary condition is satisfied^^ by placing 
for each physical vortex at location x a single image vortex with opposite circulation at 
location xR^ /\x\'^. This leads to an invariant Hamiltonian H{xi, . . . ,a;jv) as 

N,N , JV 



^ = E r^r, in(4) + ^ E HR' - r^) (n.2) 

N,N 

- ^ rj, ln(i?4 - 2R'x. ■ X, + rlT% 



47r ^ ' ' ' 47r 

N,N 

Ait 



i>j 

yf, rfj — {xi — Xj)"^ + (yi — j/j)^, and the asymmetry in the logarithmic 
terms arises because the image vortices do not move with the implied physical velocity at 
their location. In other words, the circular cylinder wall cannot be removed. The double 
sums run over all pairs i > j, i.e. there are N{N — l)/2 individual terms in these. As 
< Vi < R, the phase space is clearly bounded and its volume is {ttR^)^ . In the generic 
case of Fi with different signs we clearly have H e (— oo, +oo) due to various combinations 
of terms with — > i? or r.y — s- 0. Indeed, even at fixed H = E, say, individual terms in 



(11.2) can go to ±oo whilst adding up to a finite number. 



The first sum in (11.2) involves the usual free-space interaction term, which goes to 
-|-oo as rij — !■ for same-signed vortex pairs, and vice versa for oppositely-signed dipoles. 
As is well known, this symmetric appearance masks a quite distinct dynamical behaviour 
in these two cases, with same-signed vortices orbiting each other whilst oppositely-signed 
ones propagating along a straight line. The second sum involves a self-interaction term 
for each vortex, which leads to counter-clockwise propagation at fixed for > and 
vice versa for F^ < 0. Unlike the pair interaction terms, this term always goes to — oo as 
ri — > R, which is linked to the ever-closer approach of the vortex to its oppositely-signed 
image in this limit. In other words, the cylinder wall is a location of infinite negative 
energy for each vortex. Finally, the third sum in (UL3) involves the interaction of each 
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vortex with the images of all other vortices. Its terms become singular only if r^, rj —^ R 
and rij 0. 

One can note in passing that the dynamically active nature of the cylinder wall that is 
expressed by the self-interaction terms distinguishes the present case from the previously 
studied doubly periodic case^^'-'^^. These self-interaction terms add advection parallel to 
the wall to the dynamics, which somewhat enhances the mobility of individual vortices. 
Also, one can note that in the cylinder case propagating vortex dipolcs split up when they 
approach the cylinder wall and then propagate along the wall into opposite directions. 
Unless they collide with other vortices beforehand, the vortices would then rejoin on the 
other side of the cylinder, and again enter the interior as a dipole. Such vortex behaviour 
might in fact be relevant for vortex dynamics on beaches, where the vicinity of the 
shoreline has a similar effect as the cylinder wall^^. 



The infinities of the various terms in (IL2) occur on a set of measure zero in phase 



space volume, but they nevertheless have an impact on direct numerical simulations as 
well as on statistical theories, especially those with non-uniform phase-space measures 
(e.g. ( |ll.3| ) below). Indeed, the possibility of negatively infinite self- interaction energy is 
important even in the simplest case = F = const, although this seems to have been 
overlooked at times. For instance, Caglioti et al.^ consider the conditions for existence 
of the usual canonical partition function Z defined by the total phase space integral 



J exp(-/3i?)da;^, (n.3) 



where /3 is the usual parameter inversely proportional to the (statistical) temperature. 
They state that Z < +00 exists for fixed N and F if and only if 

/-87r 



where the finite negative range is needed in order to bound the importance of same- 
signed vortex collisions with their infinite positive energies. However, this miscalculates 
the importance of the negative infinite energies as vortices approach the wall. Indeed, Z 
for a single vortex in a cylinder is easily evaluated as 

Z=-^R'('-^'^ only if = ^ < 1, (II.5) 

1 — p* 47r 



otherwise Z < +00 does not exist. It seems that this implies that ([1.4) in Caglioti et al. 
should be replaced by the slightly more symmetrical 

which exhibits a finite temperature range also for /3 > 0. (In the particular asymptotic 



limit subsequently studied in that paper®, the re-scaled upper limit in (II. 6) still tends 
to -l-oo, so the subsequent results may well remain intact.) 

Somewhat surprisingly, this means that a point vortex system with F^ = F and bounded 
by a solid wall, if coupled to an infinite energy reservoir at positive temperature T, will 
collapse to the wall as T drops towards F^/47r, where Boltzmann's constant has been set 
to unity. Some evidence for such behaviour is given in §0 below. Of course, the point 
vortex model will lose its physical significance when this happens, i.e. the finite core 
and finite self-energy of physical vortices will become important in this limit. Also, any 



coupling to a finite energy reservoir will arrest the collapse (cf. §111). 



Now, in the remainder of this paper the following set-up will be studied. The total 
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number of vortices N = 100 will be split into Na = 4 strong vortices and Nb = 96 weak 
vortices with respective circulations 

Ta ^ ilOTT and Tb = ±2tt. (11.7) 

The net circulation in either group is zero. This is a natural constraint for physical 
situations that have arisen from localized vortex forcing, which always produces vortex 
dipoles with zero net circulation^^. These particular values have been chosen in order to 
focus on the most complex scenario, as follows. 



Because absolute energy values in (II.2) are meaningless (unlike in most classical sys- 



tems), the relevant measure of importance of the individual sums in (II. 2) comes from 



considering their variance over the phase space. Neglecting Ta^b interaction terms be- 



tween strong and weak vortices, it turns out (cf. §111 below) that the total variance of 
the terms in the double sums scales approximately as 0{N\T\), and correspondingly 
so for the weak vortices. Therefore the parameters have been arranged such that these 
variances are approximately equal. The total variance of the r\ self-interaction terms in 
(II. 2) scales as 0{NaT\) and it turns out that for the chosen set-up the numerical pre- 



factor makes this term significant in size compared with the term 0{Nj^T\). However, for 
the weak vortices the self-interaction variance 0{NBTg) is small compared to the term 
0{N^r'^). In summary, strong and weak vortices are expected to interact vigorously and 
the dynamics of the strong vortices is in addition significantly affected by the presence 
of the wall. 

It remains to consider a second invariant that arises due to the azimuthal symmetry 



of the Hamiltonian H in (11.2), namely the invariant angular momentum 

N 

2tt 



1 ^ 

^-^T.^^^- ("-8) 



i=l 



There are no other invariants, so vortex motions with N > 2 are presumably non- 
integrable in the cylinder. Unlike H, the invariance of M is not robust in the sense that a 
small disturbance of the problem (say, perturbing the cylinder wall to be elliptical) would 
destroy the invariance of M. Nevertheless, in the present case M clearly plays a role and 
needs to be considered formally on the same footing as H. It is straightforward to show 
that for the chosen set-up the strong and weak vortex contributions to the variance of 
Af are again roughly equal. 

Finally, it is noteworthy that in a set-up in which all the Ti are sign-definite, the 
conservation of M implies that the accessible phase space is bounded even without a 
cylinder wall, and this has been used to study negative temperature states of such a 



set-up using only the first term in (|IL2|). However, in the present set-up with Tj of e ither 

was 



sign this cannot be done (contrary to assertions sometimes made^'^*, where (II 
misquoted with Ti replaced by Ff 



III. A SOLVABLE TOY MODEL 

Here a toy model with exactly solvable statistical mechanics is discussed in order to 
prepare the ground for the statistical mechanics of the vortex system. The toy model is the 
one-dimensional Ising model without external magnetic field, which can be thought of as 
an assembly of independent switches. Let there be M switches and let each switch either 
be in an "up" or "down" position, with corresponding energy values Ei — ie^, where 
i = 1, 2, . . . , M and the constants €i describe the individual strength of the switches. The 
finite-sized discrete state space of the system is formed by the 2^^ different states of all 
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the switches. The total energy is 



M M M 



E^^Ei with range - ^ |ei| < £^ < + ^ |e.|. (III.l) 



4=1 i = l 



The energy extrema correspond to exactly one state each, and it is clear that there are 
only a few states with energies near these extrema. This scarcity of states is of crucial 
importance for the statistical mechanics of this system. Now, in an analogy with vortex 
dynamics each of the M switches corresponds to one of the 0{N'^) interaction terms in 



the first sum in (11.2), the other sums being disregarded. Each then corresponds to a 
pairing TiTj, and the up/down switch positions correspond crudely to the variability of 
the logarithms. Surprisingly, it will turn out that this crude analogy can already explain 
a number of features of the vortex system. 



III.l. Canonical statistical mechanics 

We first imagine the system in contact with an infinite energy reservoir at temperature T 
and inverse temperature P — 1/T. The individual switch statistics are then independent 
from each other and hence it suffices to consider the ith switch in isolation. (The anal- 
ogous statement is of course not true in the vortex system, which is a coupled A^-body 
problem.) The probabilities for the switch to be in the up/down position are then equal 
to exp(=p/3ei)/Zi respectively, where the normalization constant is the partition function 
Zi = 2cosh(/3ei). The resultant average switch energy Ui is 

C/, = -e,tanh(/3e,), (III.2) 

which shows that there is no energy equipartition between the switches unless (3 = 0, 
and that regardless of the sign of a positive temperature corresponds to negative Ui 
and vice versa. The behaviour of the switch as a function of f3ei is easily characterized: if 
P\ei \ :s> 1 (i.e. T 0+) we have Ui « — \ei\ and the switches are increasingly locked into 
their low-energy positions. On the other hand, if f}\ei\ 0+ (i.e. T — > oo) the switches 
become increasingly disordered and Ui goes to zero. A corresponding scenario unfolds 
for /9 < 0, which shows that positive and negative temperature statistics are perfectly 
symmetric here. All this applies qualitatively to the vortex system as well, though the 
near-symmetry between positive and temperature states is lost there when considering a 
mean-field theory*^^. 

The switch fluctuations can be analyzed by considering the number Ui/ti = — tanh(/3ei), 
which is the average switch position. (Identical conclusions are reached by considering 
the variance of the switch energy, which is equal to —dUi/d(3 = ef / cosh^ (/Se^ ) . ) Absolute 
values of this number near unity mark ordered, predictable behaviour, whereas absolute 
values near zero mark disorder and randomness. For uniform order increases as 
increases. Consider now the case of there being two distinct types of switches: strong 
and weak, respectively, such that [caI > |e_B|, in obvious analogy to the vortex set-up. 
Clearly, for the same value of j3 the behaviour of the strong switches is going to be more 
ordered than that of the weak switches, i.e. {Ua/^a] > \Ub/^b\: with equality holding 
only when /? = 0. Indeed, if le^l ^ les] then the strong switches can exhibit significantly 
more order than the weak switches, which illustrates the first part of Onsager's remark. 

This point about order/disorder can be made in another way by looking at the con- 
tributions to the entropy of the system that are made by the strong and weak switches, 
respectively. The usual expression for the entropy in the canonical formalism is 



Si = \i\Zi + (3Ur = ln(2cosh(^e,)) - /3ei tanh(/3ei). 



(III.3) 
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The switch entropy Si{Pei) is hence a positive even function of its argument with max- 
imum Si{0) = ln(2) at the origin, and with monotone decrease to zero with increasing 
|/3ei|. (At /3 = and /? ±00 the value of Si is compatible with the microcanonical 
definition of entropy as the logarithm of the number of permissible states.) Again, this 
means that at fixed (3 we have Sa < Sb, i-e. the strong switches contribute less to the 
entropy than the weak switches. 

A dynamical interpretation of this remarkably composite, inhomogeneous entropy dis- 
tribution is suggested by the alternative variational formulation of canonical statistical 
mechanics in terms of maximum (information) entropy at fixed mean energy. In this 
view, the strong switches tend to become more ordered, and hence contribute less to 
the entropy, because in doing so they absorb the right amount of energy to allow the 
weak switches to be as disordered as possible. This peculiar sharing-out of the energy 
leads a distribution of maximum total entropy, which then emerges as a truly composite, 
interactive feature of the system. This illustrates the second part of Onsager's remark. 



We now turn to consider the microcanonical statistics of the toy model, in which the 
total energy is fixed at a certain value E such that only states with this energy value 
are permitted and all such states are then deemed equally likely. This is the proper 
setting in which to analyze numerical simulations at conserved energy under an ergodic 
approximation, and the formalism developed here will be directly relevant to the vortex 
dynamics discussed later. 

The total energy constraint now couples the individual switch statistics. Specifically, a 
subset of switches now behaves like a subsystem in contact with a finite energy reservoir 
formed by the other switches, and the key question is to analyze the statistical mechanics 
of such a subsystem. To this end we split the M switches into a subsystem A and a 
"reservoir" B such that Ma + Mb = M and Mb '> 1, assuming from now on that 
M 1 to begin with. We have 



for all permissible states. The probability of a particular subsystem state with energy 
Ea is then proportional to the number of reservoir states with Eb = E — Ea- Now, the 
density of reservoir states per unit energy interval is that of a sum of Mb ^ 1 independent 
zero-mean random functions with finite variances, which by the central limit theorem can 
be approximated by the continuous function 



is the energy variance of the reservoir. The subscript denotes that this probability 
corresponds to a uniform distribution over all possible states. For a particular subsystem 
state with Ea we therefore obtain the probability 



III. 2. Microcanonical statistical mechanics 



E^Ea + Eb 



(III.4) 




where 




(III.5) 



Prob {subsystem state with Ea} oc po{Eb = E — Ea) oc exp 





(III.6) 



where a factor independent of Ea has been absorbed into the proportionality constant 
determined by normalization. This means that for small subsystem energies E^ <C 2a% 
the relative probabilities are described by canonical statistics with an inverse temperature 
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(3 — —E/a'^g. We note in passing that if \Ea\ <C \Eb\ then (3 agrees with the microcanon- 
ical definition of inverse temperature of the reservoir as the derivative of \tipq{Eb) with 
respect to Eb- 

Now, for larger subsystem energies the finite size of the reservoir is feh via the second 
term in (III.C), which decreases the hkehhood of such states and marks the departure from 
canonical statistical mechanics for the subsytem. (Coincidentally, the statistics of a single 
switch in contact with the finite reservoir are precisely canonical, due to the symmetry of 
its two energy levels.) Also, one can note that the probability distribution (HI. 6) admits 
a somewhat unusual maximum-entropy formulation in which the subsystem information 
entropy is extremalized subject to two constraints, namely that of fixed average energy as 
well as fixed energy variance. It is straightforward to show that such a procedure leads to 
a probability (x exp(— —jEV) for a subsystem state with energy Ea, in accordance 
with the structure found in ( [II. q ). 

The above shows how microcanonical statistical mechanics can be used to calculate 
the statistics of a subsystem in contact with a finite reservoir, even beyond the usual 
asymptotic limit Mb — > oo, E /Mb = 0(1) in which canonical subsystem statistics would 
emerge. The important reservoir energy variance cr^ in (111.5) enters the definition of (3 
and it also demarcates the finite size of the reservoir as felt by the subsystem. In analogy 
with the toy model and ( III.5| ), the reservoir variance in the vortex system is proportional 
to A'^^r^, as was previously asserted. 



IV. NUMERICAL SIMULATIONS 

Direct numerical simulations of the vortex system set-up from §H are described and 
analyzed. The comparison with statistical mechanics predictions follows in ^V| . 

IV.l. Numerical model details 



The model integrates the dynamical equations derived from (11. 1) and (II. 2), i.e. 

r,/27r 



N 



N 

E 



T'/2n 



(IV.l) 



{y'j - Vi 



As noted in 



each vortex has an image vortex with parameters 



-J' 



(IV.2) 



to satisfy the wall boundary condition. Notably, the second sum in ( IV.l ) includes a self- 
interaction term with j = i. Strong and weak vortices were chosen as in (11.7) and for 
definiteness the cylinder radius R — 5. By re-scaling (IV.l) the present simulations can be 
mapped onto simulations with arbitrary finite R and Tb provided the ratio Ir^/Fsl = 5 
remains the same. 

The numerical model itself is a standard Runge-Kutta scheme of 4/5th order with 
adaptive timestep refinement. The integrations are performed with double precision and 
in cycles of 6t = 0.01, using adaptive timestep refinement until the error is less than 
the tolerance set to 10~^. No regularization of the equations for numerical purposes 
was needed, e.g. there was no near-field smoothing to prevent large velocities, and the 
unusually low error tolerance is required to integrate safely through intermittent episodes 
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Figure 1. (a) Total energy pdf (based on 10^ random samples with uniform distribution) with 
indications of the total energy levels for the Low, Neutral, and High runs. The average total 
energy is 222 with variance 284^. (b) Joint total energy-angular momentum pdf, after some 
smoothing. The correlation coefficient between E and AI is 10~^. 

in which vortices are getting very close to each other, or to the wall. These intermittent 
episodes lead to a model performance that can vary by a factor of ten over a run. The 
total energy and angular momentum are conserved with very high precision and all runs 
were performed on a single-processor workstation. 

IV. 2. Description of three model runs 

Under the ergodic approximation the statistics of a run are determined by the invariant 
values of its Hamiltonian H — E and its angular momentum M = M, say. All runs had 
values of M very close to zero. This is because no exploration of the role of M was intended 
because of the non-robust nature of this invariant noted in §||. Three different energy 
values (denoted Low, Neutral, and High) were chosen to give runs broadly corresponding 
to regimes of positive, zero, and negative microcanonical temperature. The specific values 
for the three runs were 

= {-197,221,628}, M = {2.1, 4.1, 2.3} (IV.3) 

and they were determined as foUows^*^. A random population of 10^ vortex configurations 
was generated in which each of the hundred vortices was placed independently with 
uniform probability anywhere inside the cylinder. The H values were then computed 
from ( p^I.2| ) (this being the computationally expensive step) and a histogram was formed 
to give the probabiHty density function (pdf) for the total energy po{E), as plotted in 
figure |a. The average energy is non-zero because of the sign-definite effect of the self- 



interaction terms in (II. 2) and also because there are 0{N) more oppositely-signed terms 
than same-signed terms in the double sums. The corresponding pdf for the total angular 
momentum (not plotted) is close to a zero-mean normal distribution with variance equal 
to r|i?4/487r2 summed over all vortices. Also plotted in figure ^ is the joint pdf of total 
energy and total angular momentum pq(E,M), which indicates approximate statistical 
independence of E and M. 

At the begin of each run the four strong vortices were placed in the same positions, i.e. 




Figure 2. Inital states for the three runs with Low, Neutral, and High energy levels, from left 
to right. Large circles show the strong vortices, and black or white colour indicates positive or 
negative circulation. 



these vortices were symmetrically spaced in azimuthal angle at a common radius r — 3, 
with alternating circulation from vortex to vortex. This is very close to an exact steady 
state of four vortices in a cylinder, which occurs at radius r = (\/T7 — 4)^/'*i? « 2.96, 
and hence the strong vortices are essentially set into motion by the weak vortices. The 
positions of the weak vortices were taken from uniformly distributed random samples 
that were generated until a sample with suitable E and M was found. The resulting 
initial states are displayed in figure ^j. It is probably fair to say that it is impossible to 
guess by inspection of figure || whether and in what way the long-term statistics of the 
strong vortices will differ in these runs. 

All runs were now integrated up to t = 750, which corresponds to many hundred 
cylinder traversals of each vortex, and 3000 instantaneous states of the evolution were 
stored at time intervals of At = 0.25. Very different behaviour of the strong vortices 
was observed. The low-energy run showed a strong tendency of vortices sticking close 
to the wall, or forming short-lived vortex dipoles. The high-energy run showed a strong 
tendency to form same-signed vortex pairs that persisted a comparatively long time. 
The neutral-energy run showed a mixture of both behaviours. In all cases there were 
fluctuations around this behaviour, but the self-organization of the vortices into these 
typical patterns was conspicuously clear. 

Several quantitative diagnostics for the strong vortices were computed: their energy 
Ea (i.e. ([L2) evaluated using only the strong vortices), the distance between same- 
signed vortex pairs, and the distance between oppositely-signed vortex dipoles. The mean 
values and standard deviations of these quantities over the duration of the simulations 
are summarized in table |l|, with pdfs presented in the next section. These numbers make 
clear that significant statistical differences are indeed observed between the runs. With 
increasing overall energy E the average of the subsystem energy Ea (which initially has 
equal value in all runs) increases and the average distance between same-signed vortex 
pairs decreases. Also, at the neutral energy level (with zero microcanonical temperature) 
there is no preference between pairing of same-signed or oppositely-signed vortices. The 
average distance between oppositely-signed vortex dipoles remains roughly constant as 
E is increased and close to a value corresponding to uniformly random placement of the 
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Low Neutral High 



Energy Ea -23 ± 126 141 ± 111 269 ± 136 

predicted -64 ± 136 129 ± 117 242 ± 152 

Same-sign distance r,j 5.1 ± 2.1 4.5 ± 1.9 3.7 ± 2.0 

predicted 5.1 ± 2.1 4.5 ± 1.9 3.9± 2.1 

Opp.-sign distance r,j 4.5 ± 2.5 4.4 ± 2.1 4.6 ± 1.8 

predicted 4.6 ± 2.5 4.5 ± 2.1 4.5 ± 1.9 

Radius n 3.6 ± 1.2 3.3 ± 1.1 3.1 ± 1.0 

predicted 3.6 ± 1.2 3.3 ± 1.1 3.1± 1.1 

Table 1. Mean values and standard deviations over the duration of the simulation vs. theoretical 
predictions for various quantities. For the distances rij and ri these quantities are averaged over 
all relevant combinations of i,j £ {1, 2, 3, 4} if ri_2 ~ — IOtt and = -I-IOtt. 



vortices, vifhich is 4.5 (with standard deviation 2.1). However, the standard deviation of 
this quantity reduces with E. 

Table |l| also includes predictions based on the theory described in §0 below. These 
predictions are generally quite accurate, except for Ea in the low-energy case and for 
the same-signed in the high-energy case. Sample autocorrelation functions based on 
the numerical time series feeding into these averages were computed in order to estimate 
confidence intervals. These suggested long-lived oscillations in the low-energy E^ with 
periods of about t ^ 50 as well as a slowly decaying autocorrelation in the high-energy , 
which presumably is linked to long-lived vortex pairs in that case. Confidence intervals 
based on these estimates put the observed discrepancies at the borderline of statistical 
significance. 



V. COMPARISON WITH STATISTICAL MECHANICS 

The numerical results are compared with pdfs estimated based on microcanonical 
statistical mechanics for the whole system. The general estimation procedure is described 
in some detail, not least because it contains an important ad hoc approximation that 
significantly reduces computational cost. 

V.l. Estimation of pdfs 

The estimation of pdfs is based on phase space densities p{xi, . . . , a;jv), which for conve- 
nience in the present case are normalized such that 



pdcc^ 1, (V.l) 

where the integral is extended over the whole phase space. Relative to a chosen p the pdf 
of any phase space function $(a;i, . . . ,xn) taking real values is defined as 

The scaling properties of the Dirac-J function succinctly capture the thickness of the 
layers $ G [(j>, 0-f d0] measured in the surface integral on the right. Multiple pdfs p{(t)i, (j>2) 



are defined analogously using products of ^-functions. Pdfs defined by (V.2) can be 
estimated numerically by forming histograms of ^(x) based on random samples of a;. In 
theory, optimal convergence of such a procedure requires importance sampling, in which p 
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is used as the density of the random sample. In practice, using a uniform density coupled 
with histogram increments proportional to p is much cheaper in the present case of a 
finite phase space. The special case of a uniform density po — {ttR'^)^^ is particularly 
important, for example po{E) in figure Q has been estimated from 

PoiE) = (ttR^)-^ J 5{H - E) Ax^. (V.3) 

All pdfs calculated from the uniform density are denoted by po(')- 

The usual microcanonical density based on energy _ff e [E,E + dE] is defined as 

p,^ ^(^-^) (V4) 

and corresponding pdfs will be denoted by Pe{')- For simplicity, consideration of the 
angular momentum invariant is deferred until later. From ( |V.2| ) one obtains 

(X / (5($ - cj))6{H - E) dx^ oc poiq^, E), (V.5) 



up to an overall normalization factor. This means that pEifj^) can in principle be evaluated 
from a joint pdf po{(t>,E) based on the uniform distribution. However, for a particular 
value of E this is computationally very expensive, as most samples have to be discarded. 
On the other hand, if $ depends only on a subset of the variables then its pdf can be 
much simplified, as follows. 

Specifically, let xa and x b denote the coordinates of all the strong and weak vortices, 
respectively. Then dx^da;^ = dx^ and we consider only functions $(a;^) from now on. 
From (|V.5|) one obtains 



Pe{4>) oc / (5($ 



6{H - E) dxB 



dXAOi J ~ (j))pE{xA)dXA, 



(V.6) 



where the integral in squared brackets is the marginal density pe{xa)- Now, if it were 
the case that H = Ha(xa) + Hb{xb) (as was true in the toy model) then 

Pe{xa) oc j 5{Hb +Ha-E) dxB oc po{Eb = E - Ha), where (V.7) 

PoIEb) oc J 5{Hb - Eb) dxB (V.8) 
is the pdf of the weak vortex energy Hb in the uniform distribution. The last term in 



(V.7) means that the function po{Eb) should be evaluated at Eb — E — Ha- Clearly, 
all states xa with the same energy Ha are now equally likely. Because pe{xa) depends 
only on E — Ha this is a huge simplification, as poiEB) can be computed once and for 
all from the uniform distribution. 



However, the Hamiltonian (II. 2) does not fit into this category, as we have 



H{xa,xb) = Ha{xa) + Hb{xb) + Hi{xa,xb) (V.9) 

where the 'interfacial' energy Hj consists of terms involving both strong and weak vortex 
circulations. In principle, this means that pe{xa) does not depend solely on E — Ha- 
Calculating pe{xa) for all xa directly would be very expensive, as even in the present 
case this would require a look-up table in eight dimensions. Instead, a much simpler ad 
hoc approximation for pe{xa) as a function oi E — Ha is used, which for Hi = reduces 



to (V.7). The approximation is 

Pe{xa)ocpo{Er^E-Ha{xa)), (V.IO) 
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where 



Po{Er) (X J j 5{Hr- En)dxA<^xi 
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is the pdf of th e 'reservoir' energy Hr = Hb + Hi based on the uniform distribution. 
Unhke in ( V.8 ), the double integral is necessary here because Hr depends on both xa 



and xb- This approximation effectively assigns the same probabilitiy to all states xa 
that have the same strong vortex energy Ha{xa)- Indeed, the assigned probability can 
be shown to be the average probability over all states with the same Ha- 
Based on this approximation, the pdf for any <^{xa) is now given by 



Pe{4>) « y <5($ - 4>)pq{Er = E-Ha) Axa 



The pdf of Ha in particular simplifies further to 

Pe{Ea) « po{Ea)po{Er = E- Ea). 



(V.12) 



(V.13) 



The additional consideration of the second (non-robust) angular momentum invariant 
(Q.S) is straightforward, especially as M = Ma + Mb holds exactly. The microcanonical 
density becomes 

S{H - E)S{M - M) 



Rem 



J 5{H - E)5{M - M) dx^ '' 



the marginal density 

Pem{xa) oc \ ^ii- E)5{M - M) Axb. 
and the joint 'reservoir' pdf based on the uniform distribution 



Pq{Er, Mb) oc / / 5{Hr - Er)5{Mb - Mb) Axa^xb 



(V.14) 



(V.15) 



(V.16) 



The same approximation procedure for ( V.15| ) as before then gives 



PEM{(j}) OC / (5($ - (I))po{Er = E- Ha, Mb = M - Ma) Axa 



(V.17) 



The predictions of the statistical mechanics theory are hence pdfs estimated based on 
(V.16) and (V.17). For comparison with figure |l| the functions pq{Er) and pq{Er, Mb) 
are plotted in figure |^. 

Despite the somewhat opulent appearance of ( |V.16 ) and ( V.17 ) , the practical procedure 
for estimating pdfs is actually disarmingly simple. A sample of 10^ states was generated 
using the uniform distribution and the corresponding values of H, Ha, Hr, M, Ma, Mb 
were computed and stored in lists, as were the coordinates xa of the four strong vortices. 
Computing H is by far the most expensive step here. Histograms based on these lists were 
then used to estimate (V.16). For any function ^{xa) to be investigated a corresponding 
list of values was then computed from the stored coordinates. This list together with a 
look-up table for the histogram increments oc po{Er, Mr) at the shifted arguments was 
then used to estimate (V.17) at fixed E and M. It is worth stressing that only a single 
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Figure 3. (a) Pdf of the reservoir energy Eu = E — Ea, where Ea is the energy of the 
strong vortices, based on 10^ random samples with uniform distribution. The average reservoir 
energy is 109 with variance 229^. (b) Joint reservoir energy-angular momentum pdf, after some 
smoothing, where Mb = M — Ma and Ma is the angular momentum of the strong vortices. 
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Figure 4. Pdfs of Ea, the energy of the strong vortices, for the three runs. Thin lines are 
simulation results, thick lines are theoretical predictions, and the squares denote po{Ea)- 

large sample based on the uniform distribution is needed to describe the statistics of the 
system at arbitrary E and M. 



V.2. Comparison with model results 

The pdfs for Ej^ (i.e. the energy of the strong vortices) are plotted in figure |[ Throughout, 
thin lines denote pdfs estimated from histograms taken from the direct numerical simu- 
lations and thick lines denote theoretical predictions based on (V.17). Also shown in the 
middle panel is po{Ea), which corresponds to random placing of the vortices with uniform 
distribution. As could be expected, in the neutral energy case this gives a reasonable first 
approximation, though still a less accurate one than pem{Ea)- The mean energy of the 
strong vortices increases as E increases and this can clearly be predicted quantitatively 
from the theoretical predictions. In most cases, there is very good agreement between the 
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Figure 5. Pdfs of rij between strong vortices of the same sign, for the three runs. Thin lines are 
simulation results, thick lines are theoretical predictions, and the squares denote po{rij). The 
simulation results are averages over the pdfs of ri2 and r34. 
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Figure 6. Pdfs of Vij between strong vortices of opposite sign, for the three runs. Thin lines are 
simulation results, thick lines are theoretical predictions, and the squares denote po{rij). The 
simulation results are averages over the pdfs of ri3,ri4,r23, and r24. 



theoretical and simulation statistics, not only in terms of accurate prediction of low-order 
moments, but also in the prediction of the non-Gaussian shape of the pdfs. 

Figures |5|-^ show the pdfs for the distance between same- and opposite-signed strong 
vortices, respectively. The tendency for same-signed vortices to cluster at smaller 
with increasing E is evident in figure ^. On the other hand, the indifference to E of the 
average between opposite-signed vortices that was noted in table |^ masks notable 
changes in the pdf that are evident in figure ^. These are well captured by the theoretical 
predictions. 

Finally, figure ^ shows statistics for r^, the vortex distance from the origin. This quan- 
tity is interesting due to the influence of the self-interaction terms in ([1.2), which, as 
noted in §||, significantly affect the dynamics of the strong vortices. For neutral and high 
energies the pdf of settles down to a shape quite close to the uniform distribution 
(except near the wall — R), which is indicated by the squares in the middle panel. 
However, the first panel in figure ^ shows how the strong vortices tend to accumulate 
at the cylinder wall for very low energies. This is essentially a similar statistical effect 
as the formation of opposite-signed vortex dipoles in the first panel of figure ^, but the 
wall effect is clearly more pronounced in its pdf. Interestingly, the effective temperature 
estimated as d\npo{E)/dE from figure |l| at E — —197 gives « 0.01 for this run. The the- 
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Figure 7. Pdfs of Vi, the distance from the origin of the strong vortices, for the three runs. 
Thin hnes are simulation results, thick hues are theoretical predictions, and the squares denote 
Po{ri). The simulation results are averages over the pdfs of ri,r2,r3, and r4. 



oretical upper limit in (II.6) for the possible existence of canonical subsystem statistics 
in the present case gives « 0.013, which confirms that this run is close to a collapse to 
the wall. Finite point vortex statistics in a collapse case would rely entirely on the finite 
size of the reservoir formed by the small vortices. 



VI. CONCLUDING REMARKS 

The theoretical predictions based on an ergodic approximation for the whole vortex 
system were seen to predict the pdfs of many, though not all, descriptive variables of the 
strong vortex subsystem with a surprisingly high degree of accuracy. Some simulation 
averages converged only very slowly and significantly longer integrations could test the 
quality of the ergodic approximation in these cases. On the other h and, non-ergodic 
observations might also be due to the crucial approximation leading to ( V.17 ), which was 
necessary to estimate the theoretical pdfs at affordable computational cost. Otherwise, 
the prediction procedure was remarkably simple and cheap, relying only on a single 
random sample of vortex configurations to predict pdfs for all values of total energy E 
and angular momentum M. 

The toy model, the direct numerical vortex simulations, and the theoretical predic- 
tions all corroborated Onsager's crucial insight that strong vortices will exhibit amplified 
statistical tendencies compared to weak vortices, and hence will be more predictable in 
a negative temperature state. It is the strongly unequal circulation strengths that allows 
the flow to organise itself in this inhomogeneous manner. It can be noted that in terms of 
entropy as a measure of accessible phase space volume, clustered vortices always present 
a low-entropy state. The crucial point is that in a negative temperature state the low en- 
tropy of the clustering strong vortices is more than compensated for by the high entropy 
of the freely roaming weak vortices. 

It is intriguing to note that on the level of individual vortex dynamics there is a 
smooth transition from positive to negative temperature behaviour, exemplified by the 
smooth po{E) in figure |^. By contrast, in a coarse-grained picture two vortices with 
opposite circulations that are close together cancel each other out and hence disappear 
from view. This illustrates why solutions to mean-field theories (such as the Sinh~Poisson 
equation^) have a characteristic cut-on behaviour as /3 < 0, because only same-signed 
vortex clustering is observable in the coarse-grained variables of these theories. It is also 
noteworthy that the near-collapse of the vortices to the wall in the present low-energy 
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case (which is linked to the theoretical upper bound in ( p^I.GD ) docs not occur in the 
Sinh-Poisson equation because there the scaling has been arranged ab initio to render 
the wall-induced self-interaction energies negligible. 

The present set-up seems to have been the most complicated to study, i.e. the strong 
vortices interact vigorously with both the wall and with the small vortices, and the 
small vortices themselves form only a finite energy and angular momentum reservoir. As 
noted before, the latter point puts the strong vortex statistics beyond the reach of the 
usual canonical theories. In other words, whilst the overall energy regime can be broadly 
classified by the sign of the usual statistical mechanics temperature, the temperature 
concept alone is not sufficient to make quantitative predictions. The modified maximum 
entropy principle suggested by the Gaussian pdfs of the toy model (in which entropy 
was maximized subject to both a fixed mean energy and a fixed energy variance) might 
allow analytical progress to be made here. Another possible avenue for future analytical 
progress is an asymptotic exploitation of the small parameter Irs/Fyil. 

It is tempting to generalize the present results (which directly apply only to well- 
separated finite-size vortices) to continuous vorticity distributions by letting iV — s- oo 
in some way. As is now well known, the relevant scaling must keep NT — 0(1) in this 
limit. However, this implies that the microscopic vortex mobility due to the induced 
velocity by the nearest vortices at average distance oc R/y/N would decrease to zero as 
V^/N oc 1/y/N. Therefore the microscopic vortex system equilibration time (which is 
implicitly assumed here to be small compared to the observation period) goes to infinity 
in this limit. This problem is directly observable, say, in the simplest case of many 
Ti = const, vortices in a cylinder, for which exact solutions of the N ^ oo limit are 
known*. Simulations starting from non-equilibrium initial conditions clearly show that 
these solutions are practically unattainable due the lack of vortex mobility. On the other 
hand, there is evidence'^ that the slow evolution of well-mixing large-scale fiows can be 
approximated by statistical point vortex theory. Notably, random forcing also helps in 
this context, as it too increases the mobility of the vortex population. This suggests that 
a limit in which Na = const, and Nb — > oo could perhaps be used to model the behaviour 
of strong vortices surround by a sea of "filamentary vortex debris" , because the stirring 
by the strong vortices could provide the essential mobility for the debris. 

Finally, the prediction tools developed in this paper could be used to study an interest- 
ing 'inverse' problem: from observing only the strong vortices, can one deduce the number 
and strengths of the unobserved weak vortices? This would provide a nonlinear method 
for estimating unobservable sub-gridscale data, perhaps with applications in geophysical 
fluid dynamics. 
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